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Coupled nonlinear oscillators can exhibit a wide variety of patterns. We study the Brusselator as 
a prototypical autocatalytic reaction-diffusion model. Working in the limit of strong nonlinearity 
provides a clear timescale separation that leads to a canard explosion in a single Brusselator. In 
this highly nonlinear regime, it is numerically found that rings of coupled Brusselators do not follow 
the predictions from Turing analysis. We find that the behavior can be explained using a piecewise 
linear approximation. 


I. INTRODUCTION 

The emergence of synchronized patterns in 
systems of coupled oscillators is a phenomenon 
that controls the function of a wide variety of 
systems ranging from biological clocks to opto¬ 
electronic circuits p]. Oscillators can synchro¬ 
nize in phase or out of phase depending on the 
nature of the coupling, and the inherent com¬ 
plexities of the individual oscillators. A partic¬ 
ularly intriguing question is the impact of frus¬ 
tration on phase and frequency patterns of cou¬ 
pled nonlinear oscillators [2j [3J. Recent experi¬ 
ments in a model of the classic Turing ring of 
diffusively coupled synthetic “cells” [H have ob¬ 
served a variety of synchronized patterns with 
complex phase relationships, and spatially het¬ 
erogenous patterns where the oscillation of one 
or more units was suppressed 0. The experi¬ 
mental model consisted of microfluidically pro¬ 
duced surfactant-stabilized emulsions in which 
ordered rings of aqueous droplets containing the 
Bclousov-Zhabotinsky (BZ) oscillatory chemi¬ 
cal reactants are dispersed in oil[5]. The domi¬ 
nant phase relation between oscillators was anti¬ 
phase (7T out of phase) even in odd-numbered 
rings which are geometrically frustrated [2], and 
an explosion of patterns occurred at a ring size 
of five [5]. 

Our work is motivated by these experiments 
and the broader question of the emergent be¬ 
havior of a set of coupled, strongly nonlin¬ 
ear oscillators in a geometrically frustrated ar¬ 
ray. In such an environment, coupled oscil¬ 


lators can relieve the frustration by synchro¬ 
nizing in phase, suffering oscillator death[H], 
or exhibiting complex spatio-temporal patterns 
that accommodate the frustration but preserve 
the anti-phase relations. We find that in the 
regime of strong nonlinearities, the phase space 
is dominated by complex spatio-temporal pat¬ 
terns whereas weakly nonlinear oscillators syn¬ 
chronize in phase or suffer oscillator death. 

The systems we study are Turing rings where 
the chemistry of the “cells” is modeled by one 
of the simplest autocatalytic reactions: the 
Brusselator [7]. These Brusselators are cou¬ 
pled through the diffusion of the inhibitory 
species, which mimics the diffusion in a num¬ 
ber of situations including the microfluidic 
emulsion [3], Coupled Brusselators have been 
studied extensively [5] [j7] PH], however, attention 
has been mainly focused on the regime where 
the nonlinearity is weak|2] or moderate. In the 
weakly nonlinear regime, there is a large sep¬ 
aration of time scales with the activator be¬ 
ing the fast species |9] and a canard explosion 
leads to relaxation oscillators. We focus on the 
regime of strong nonlinearity , which is also char¬ 
acterized by a large separation of time scales 
and a canard explosion|ll| leading to relaxation 
oscillations with a qualitatively different form 
of the limit cycle [Ij. In this regime, the in¬ 
hibitor is the fast species, and the slow variable 
is the total concentration of chemicals: activa¬ 
tor plus inhibitor. We show that coupling Brus¬ 
selators in this strongly non-linear regime leads 
to a wide variety of spatio-temporal patterns 
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that are characterized by a strong preference 
for neighboring oscillators to be 7r out of phase. 

We characterize the attractor space of rings 
of oscillators, coupled via the inhibitory species, 
through numerical analysis and application of a 
piecewise linear approximaiton (PLA) that is 
valid in the large nonlinearity regime. Anal¬ 
ysis of the model based on PLA shows that in 
this highly nonlinear regime the dynamics of the 
coupled Brusselators can be represented as os¬ 
cillators traveling on branches of the limit cycle 
separated by near-instantaneous jumps between 
these branches. In one branch of the limit cy¬ 
cle, the oscillators are repulsively coupled while 
on the other they are attractively coupled. The 
sequence of jumps between the branches charac¬ 
terizes the patterns observed numerically. This 
representation, in terms of jumps provides an 
understanding the rich attractor space in the 
highly-nonlinear regime, and we suggest that 
this technique will be useful for other reaction- 
diffusion systems with strong nonlinearities p~2] . 

This paper is organized as follows. First, we 
present an analytical investigation of a single 
Brusselator with emphasis on the strongly non¬ 
linear regime. We perform a change of variables 
that clearly demonstrates the separation of time 
scales, and present a review of the standard 
Turing analysis on a ring of N r Brusselators, 
using this representation. This is followed by 
a numerical investigation, which examines the 
phase and frequency synchronization properties 
of rings with N r ranging from 2 to 5. Finally we 
apply the PLA to the Brusselator in the regime 
of large nonlinearity in order to understand the 
origin of the frustration in W-rings, and the 
patterns that emerge in numerical solutions. 


II. THE SINGLE BRUSSELATOR 


In this section, we review the properties of the 
single Brusselator in order to clearly identify the 
regions of weak and strong nonlinearities and 
the time-scale separation that drives the physics 
in the regime of large nonlinearity. 

The Brusselator [13] is a model of an auto 
catalytic chemical reaction that was first pro¬ 


posed by Ilya Prigogine at the Free University 
of Brussels in the 1960s. It is a two species 
model, with an activator X and an inhibitor 
Y. There are several ways to write the set of 
reactions [T^HlF] . We choose the representa¬ 
tion below m H4] since this form is particu¬ 
larly amenable to a systematic expansion about 
the mean field limit [131 . and we would like to 
explore these effects in the near future. 


A X + A 
X -> 0 

X + B -a Y + B 
2X + Y -a 3X 


N 

n x 

bn x (1) 
cN~ 2 n 2 n y 


In this representation, the populations of the 
A and B species are constant in time, and 
their only role is to set the reaction rates[TB;. 
The first reaction describes the creation of X 
molecules with a rate proportional to N, the 
number of A molecules. The first two reactions, 
in isolation, ensure that the mean number of X 
molecules over time is given by A, which acts as 
a measure of the system size. The other param¬ 
eters appearing in the reactions are as follows: b 
is the rate of exchange of activator to inhibitor, 
c controls the rate of the nonlinear reaction in 
which the inhibitor transforms back into the ac¬ 
tivator, and n x and n y are the number of X and 
Y molecules, respectively. 

A well-mixed system, in the limit of large sys¬ 
tem size ( N ), is described by well-known rate 
equations in which the rates of change of the 
concentrations of A' (x) and Y (y) are obtained 
by using the law of mass action. The rate 
equations can also be obtained systematically 
through a Van-Kampen expansion[13|. and that 
approach demonstrates the mean field nature of 
the equations: 


x = 1 — x(l + b — cxy) 
y = x(b — cxy) (2) 

where x = X/N and y = Y/N In this paper, 
we focus exclusively on the mean field equa¬ 
tions. The presence of a rich attractor space 
and metastability suggests that intrinsic fluctu¬ 
ations, neglected in the mean field limit, can 
modify patterns qualitati velyJU, H5] . 
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The above set of nonlinear first order differ¬ 
ential equations (Eq. |2j) has a unique fixed point 
at ( x,y) = (1, (8 + c+l)/c). Where S = 6— c— 1 
is the control parameter for the Hopf instabil¬ 
ity in the system. By examining the behavior 
of these equations near this fixed point one can 
show that the parameter space of <5 and c is bro¬ 
ken up into the four regions as shown in Figure 
[l] In regions (1) and (2) where 6 < 0, the fixed 
point is stable: in regions (3) and (4), the fixed 
point is unstable. The eigenvalues of the Jaco¬ 
bian matrix are complex in regions (2) and (3) 
indicating the fixed point lies at the end of a 
stable and unstable spiral, respectively. In re¬ 
gion (3) a stable limit cycle emerges from the 
fixed point with an amplitude that grows as y/S 
for small S. 



FIG. 1: (Color Online) Regions of stability 
about the fixed point in b-c parameter space. 
In regions (1) and (2) the fixed point is stable 
while in (3) and (4) it is unstable. 
Additionally, regions (2) and (3) have 
oscillatory behavior around the fixed point. 


As c is increased beyond some threshold 
value, however, the dynamics of the Brusselator 
changes qualitatively: it enters a regime charac¬ 
terized by a large separation of time scales. This 
difference in time scales becomes apparent un¬ 


der the change of variable u = y and v = x + y. 
In this representation, the rate equations be¬ 
come. 


u = (v — w)([l + <5] + c[l — u(v — u)]) 
v = 1 — v + u (3) 


It can be easily seen from these equations that 
for c >> 1, u is the fast variable and v is the 
slow variable: x + y ~ constant when u ^ 0. In 
this regime, the Hopf bifurcation acquires a sin¬ 
gular character [IT], and at S > 1/c, the system 
passes though a canard point where the ampli¬ 
tude of the limit cycle jumps sharply. Figures 
2apb show the shape of the limit cycle for a 


fixed S as a function of c. Fig. [2b] compares am¬ 
plitudes for two different close values of S at a 
large c. The dependence of the canard explosion 
on the parameter c can be seen in Figures [2d| It 
is important to note here that in mean field the 
single Brusselator is a monostable system with 
its limiting behavior either being a fixed point, 
the small limit cycle or the large limit cycle. 


III. RING OF COUPLED 
BRUSSELATORS 


The previous section provided an overview of 
the properties of the single Brusselator in the 
regime of high non linearity. In this section, 
we extend our analysis to a ring of Brusselators 
coupled through Y. the inhibitor species. To in¬ 
vestigate this system, we introduce a ’’hopping” 
term for the inhibitor to the list of reactions, 
which translates to a diffusive coupling term in 
the rate equations. The list of reactions now 
becomes: 


A -a Xi + A 

N 

Xi —> 0 

n Xi 

Xi + B Yi + B 

bn Xi 

2X, +Yi -»■ 3Xi 

cN~ 2 n^..n yi 

Yi Yi+ 1 

dn Vi 

Yi -A Y^ i 

d n Vi 


( 4 ) 
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(a) (Color) Nullclines (dashed blue - u, dashed 
green -v), velocity field (black) and limit cycle 
(red) for c = 1, (5=1. There are two branches of 
the the u nullcline that asymptotically 
approach each other. The limit cycle grows out 
of a Hopf bifurcation of the fixed point 
(magenta dot). 



t 

(c) (Color Online)Blue (solid): Time series of u 
for c=10 <5=0.3: this is below the canard point. 
The limit cycle has a period of about 2.5 and 
an amplitude of about 0.7. Green (dashed): 
Time series of u for c=10 5=0.4: this is above 
the canard point. The limit cycle has a period 
of about 5.2 and an amplitude of about 4. 



(b) (Color) Nullclines (dashed blue - u, dashed 
green -v), velocity field (black) and limit cycle 
(red) for c = 50, 5=1. The asymptotic 
approach of the two branches of the the u 
nullcline is clearer at this value of c, and the 
limit cycle grows out of a singular Hopf 
bifurcation of the fixed point (magenta dot). 
For this value of c, the limit cycle follows the 
nullclines of u very closely. 



(d) (Color Online)Green (dash - dot), Blue 
(solid) and Red (dotted) are the scaled 
amplitudes of the limit cycle for c=l, 10 and 
100 respectively. These scaled amplitudes are 
plotted against 5 = b — c — 1, and show that 
c = 1 has no canard explosion in this range of 5 
while c = 10 and 100 do. 


FIG. 2 
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where i labels the individual Brusselators. The 
corresponding rate equations are: 

Xi = 1 - Xi{l + b - cXiyi) (5) 

Vi = Xi{b- CXiyi) + d(y i+ i - 2y t + t/*_i) 

Where d is the parameter that controls the cou¬ 
pling strength. In the u and v variables, the 
equations are: 

Ui = (vi - Ui)( 1 + < 5 ) ( 6 ) 

+ c[(Vi - Ui) - Ui(Vi - Ui) 2 } 

T d)\li -(_i 2 Ui + Ui—l) 

Vi = 1 - Vi + Ui + d(u i+ 1 - 2 Ui + Ui- 1 ) 

The original diffusive coupling is realized in 
this representation as a coupling of the fast vari¬ 
ables, u, and an additional driving term for the 
slow variable, v. 

A. Linear stability analysis of Ring 


Where now a q is the growth rate, q is the 
wavenumber and s = l,...,N r is the oscillator 
number. Substituting this form into Eq. [8j 



/ e ig(s-i) _ 2e *9s + e *90+i) 0 \ / u \ 

^ e i 9 ( s -i) _ 2e iqs + e iq(s+i) oj yvj 


The allowed values of q are q = with n 
an integer, and a q is determined by: 



with J n given by: 

f _ (l + 6-4dsm 2 (^) -(1 +S + c)\ 
Jn ~[ l-4dsin 2 (^)" -1 )‘ 


In order to understand the emergence of pat¬ 
terns in a ring of N r Brusselators we perform a 
standard linear stability analysis of the homo¬ 
geneous fixed point; that is the fixed point of 
a Brusselator without coupling. The analysis 
is summarized here for completeness. We lin¬ 
earize each Brusselator about the fixed point at 
(ui,Vi) = + i±f±^). For the uncou¬ 

pled system, the deviation from the fixed point 
obeys the equation: 


Stability of the fixed point requires that the real 
part of both eigenvalues is negative which is sat¬ 
isfied when Tr^jj is negative and det^jj is 
positive for all n. Thus the requirements for 
stability are: 

S — 4dsin 2 (^) < 0 (12) 

c — 4d(c + S) sin 2 (^) > 0 (13) 


u 

v 



(7) 


Where J is the jacobian of this system evaluated 
at the fixed point: ^ ~ (! + <> + c)\ 


Adding the diffusive coupling leads to: 


u 

V 


(J + D) 


( 8 ) 


where D is the coupling matrix for u and v We 
now look for solutions to the linear equation 
with the form: 


U = U g 


e <T q t e iqs 


(9) 


We are interested in cases where the coupling 
d « 1, and in this regime, the fixed points are 
stable for 6 < 0. For S > 0 the fastest grow¬ 
ing mode is the one for which sin 2 (^) = 0, 
therefore, the q = 0 mode. This suggests that 
the Brusselators in a ring will always oscillate 
in phase with one another. In the following sec¬ 
tions, we present results of numerical simula¬ 
tions for rings with N r = 2 — 5. We have care¬ 
fully investigated the phase-diagram in d — 5 
space for N r = 2, 5. For the N r = 2 system 
we focus on explaining the patterns using the 
piecewise linear approximation. For N r = 5 we 
emphasize the complexity of the patterns that 
form and the regions that they form in. 
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IV. NUMERICAL RESULTS 

A. Two Oscillators 


In this section, we present a numerical study 
of the synchronization of two Brusselators, fo¬ 
cusing primarily on the large c regime. From 
the linear stability analysis, the oscillators are 
always expected to synchronize in phase. Sim¬ 
ulating this system in the mean held limit we 
find that when c = 1 that two Brusselators will 
synchronize in phase with one another as the 
Turing analysis predicts. In the c> 1 regime, 
we find that when the system is above the ca¬ 
nard point - in the large limit cycle - the Brus¬ 
selators oscillate out of phase with one another, 
while below the canard explosion they oscillate 
in phase with one another. 

For N r = 2, we observe four different types 
of behavior - homogeneous stable fixed points, 
inhomogeneous stable fixed points (oscillator 
death), in phase oscillations and out of phase 
oscillations. We notice that if the coupling- 
strength d is above a critical value d c then the 
Brusselators experience oscillator death [6j: an¬ 
other, nontrivial pair of fixed points becomes 
real and stable through a pitchfork bifurcation 
for d > d c = ( e+ ^j _ 2 )2 ■ This type of behavior, 
however, is not the focus of this paper. There 
are also regions where the behavior is depen¬ 
dent on initial conditions, these are enumerated 
in Table [I] It is in these regions that hysteresis 
can be observed. A complete portrait of the at¬ 
tractor space in d — 5 space will be discussed in 
the context of the PL A, and is shown in Fig. [Toj 

Summarizing, when the system is in the small 
limit cycle, in-phase oscillations occur, as pre¬ 
dicted by the linear stability analysis. When 
the system is in the large limit cycle, the lin¬ 
ear stability analysis fails for some regions of 
parameter space, and out of phase oscillations 
can be observed. These out of phase oscillations 
are expected to lead to complex patterns when 
the oscillators are arranged in a geometry that 
is frustrated: the oscillators cannot be n out of 
phase with all of their nearest neighbors. From 
here on we only look at Brusselators with large 


Region 

Possible states 

1 

Homogeneous fixed point 

2 

Homogeneous fixed point, 
Out of phase oscillations. 

3 

Homogeneous fixed point, 

Oscillator death. 

4 

Out of phase oscillations. 

5 

Oscillator death. 

6 

Out of phase oscillations, 
In phase oscillations. 

7 

Oscillator death. 

In phase oscillations. 

8 

In phase oscillations. 


TABLE I: States observed in a numerical 
investigation of the 2 - ring. A complete phase 
diagram with numbered regions is found in 


Fig. 10 


nonlinearity: specifically, c = 50 and below the 
coupling at which they stop oscillating (d < d c ). 


B. Three Brusselators 

A triangular ring is the smallest to exhibit ge¬ 
ometric frustration: in the large limit cycle all 
Brusselators prefer to be exactly tt out of phase 
with their nearest neighbors, which is impossi¬ 
ble in this geometry. In contrast to the behavior 
of phase-coupled Kuramoto oscillators [IS], nu¬ 
merical solutions to the mean Field equations 
show that the Brusselators do not relieve their 
frustration by adjusting the phase difference to 
be 27 t/ 3. Instead, the 3—ring adjusts to the 
frustration by maximizing the number of i r out 
of phase interactions: synchronizing two of the 
Brusselators in phase and one Brusselator out 
of phase with the other two as shown in the 
schematic in Figure [3] Additionally, if the cou¬ 
pling is weak enough the Brusselators are also 
observed to oscillate in phase. In experiments 
on hexagonal arrays of oscillators [50! with the 
BZ chemistry, two patterns are observed. At 
low diffusive coupling the 27 t/ 3 pattern appears, 
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and at strong coupling a s 07 t[ 5 ] pattern is ob¬ 
served, which also maximizes the number of 7r 
out of phase oscillators. In this pattern how¬ 
ever the central droplet stops oscillating while 
the surrounding six droplets oscillate 7 r out of 
phase. 



There exist regions in 6 — d space where 
the Brusselators exhibits multistability of these 
modes. For example at S = 1, d = 0.02 modes 
1, 2 and 3 are all stable periodic solutions ac¬ 
cessible from different initial conditions. 


TT 

• • 



FIG. 3: (Color Online) Modes of oscillation in 
a triangle of 3 Brusselators. Phase differences 
between Brusselators are marked. Diagonal 
lines represent a relative phase. 


(a) Mode 1 (b) Mode 2 



0 0 0 0 




C. Four Brusselators - Multi Stability 


(c) Mode 3 (d) Mode 4 


Expanding the ring to 4 Brusselators elim¬ 
inates the geometrical frustration associated 
with odd numbers of oscillators that prefer to 
be tt out of phase. Numerical results, however, 
reveal a complex landscape of attractors with 
multiple modes of oscillations occurring at a 
fixed set of model parameters. 

Mode 1: As one would expect from the re¬ 
sults of the 2 Brusselator system there is a mode 
where all Brusselators are out of phase with 
their neighbors. 

Mode 2: The first additional mode is one 
where each Brusselator has one in phase neigh¬ 
bor and one out of phase neighbor, this leads 
to the sides of square oscillating out of phase 
with each other as can be seen in the cartoon in 


Figure 4b 


Mode 3: The next mode of oscillation ob¬ 
served is characterized by two Brusselators that 
have stopped oscillating on opposite corners 
complimented by two out of phase oscillating 
Brusselators on the remaining corners as can 
be seen inl4cl 

Mode 4: Finally there exists a mode where 
all four Brusselators oscillate in phase. 


FIG. 4: (Color Online) Modes of oscillation for 
a ring of 4 Brusselators. Diagonal lines 
represent a relative phase. 


D. Five Brusselators - Frustrated patterns 


The geometry of this odd numbered ring 
frustrates the system, as for the three-oscillator 
system since it is impossible for all the Brus¬ 
selators to be exactly tt out of phase with 
their neighbors. Unlike the 3-ring, however, 
we encounter an explosion of synchronized 
patterns for the 5-ring. The oscillating patterns 
can be characterized by five distinct relations 
as shown in Table [TT] 


Starting from different initial conditions these 
modes occur with varying probability in differ¬ 
ent parts of the S — d phase space as outlined 
in Fig. 5a- 5e Furthermore, when in mode 2 
the oscillators can exhibit different frequencies 
from each other. Specifically, the synchronized 
























Mode 

Num¬ 

ber 

Description 

Schematic 

1 

Nearest neighbors have a 
phase difference of 47r/5. 
This creates a star like fir¬ 
ing pattern 

W 

2 

One pair of nearest neigh¬ 
bors is in phase and one 
pair of second nearest 
neighbors are in phase, 
leaving one to oscilate by 
itself. 

) 

(Q 0 
(• •) 

3 

One pair of second near¬ 
est neighbors are in phase 
this traves around the 
ring. 

®\ 

a. w 

4 

One pair of second near¬ 
est neighbors are in phase 
and the the remaining 3 
are in phase. 

l^\ 

5 

All are in phase 



TABLE II: Description of phase relationships 
in observed inodes in a ring of five Brusselators 


patterns are characterized by a frequency ratio 
of the three groups of oscillators participating 
in Mode 2 (Table [II]) , but the ratios change dis- 
continuously as S is varied at a fixed d in the re¬ 


gion of the phase diagram indicated in Fig. 5a 


[5c| This behavior is reminiscent of phenomena 


such as the Devil’s staircase encountered in dis¬ 
crete systems [21]. The piecewise linear model 
discussed in the next section offers insight into 
the multitude of patterns observed in rings of 
the strong relaxation oscillators: Brusselators 
at large c. 


V. PIECEWISE LINEAR 
APPROXIMATION 

A. Approximating a single Brusselator 

In this section, we analyze the dynamics of 
W-rings at large c using a piecewise linear ap¬ 
proximation, PLA, to gain some insight into 
the multitude of out-of phase synchronized pat¬ 
terns. Specifically, we are trying to under¬ 
stand the origin of the out of phase oscillations 
and how that gives rise to the multistability of 
modes. The PLA approximation has been used 
in previous studies to analyze strong relaxation 
oscillators, such as the Van der Pol oscillator [12] 
and the Tyson-Fife model [22]. This has proven 
to be an effective method for studying relax¬ 
ation oscillators with different coupling mech¬ 
anisms including diffusive coupling [12] and de¬ 
lay coupling [23], and in systems with two[T2], 
three[?4] and four oscillators]?? . 

We begin our analysis with the single Brus¬ 
selator described in terms of the fast-slow vari¬ 
ables: it and v. 

^u=f(u,v) (14) 

v=g(u,v) (15) 

where 

/(it, v ) =(v — it) ^ —— + 1 — u(v — u)^j (16) 

g(u, v) =1 — v + u (17) 

We can see from Fig. [2b] that the large 
limit cycle consists of segments that closely 
follow the the nullclines of it, the fast vari¬ 
able, with near instantaneous jumps between 
the branches of the nullclines. Specifically, 
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(a) 


(b) 
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FIG. 5: (Color Online) Probability heat maps for the modes of oscillations of the 5-ring, (a)-(e) 
correspond to Modes 1 through 5 in Table [TT} (f) shows the probability that an initial 
configuration of identical Brusselators have a final state where they oscillate with different 
frequencies. Color indicates how often that mode of oscillation was observed in an ensemble of 50 
randomly chosen initial conditions. Blue indicates that the mode of oscillation was rarely or never 
achieved for that value of d and 5, while red indicates that the mode was achieved by a majority 
of the randomly chosen initial conditions. 
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-5 0 5 

5 


FIG. 6: (Color Online) Full phase diagram for 
the 5-ring illustrating the diversity of patterns. 
Each color represents a different combination 
of possible modes of oscillation achievable from 
different initial conditions. 


the limit cycle climbs up the line it = v un¬ 
til it falls off at (it, u) = ( u max ,u max ) where 

Umax = (c+ ig 2) • We determine u m ax from 
the maximum of the x nullcline in the x , y 
representation. After reaching u ma x, the os¬ 
cillator jumps to the left most branch of v = 
it + ,<5+ e+ c ~ and starts following that branch 
downwards at coordinates (u,v) = {\{u max — 
V u max -4(5+1 + c)/c), u max ) until it falls off 


this branch at the minimum of this nullcline at 
(it, it) = (\/(S + 1 + c)/c, 2 a /S + 1 + c/c). 



FIG. 7: (Color) Null dines of u from the full 
model (blue), and from the piecewise linear 
approximation (red) 

To capture the essence of these dynamics, 
we approximate the trajectory by restricting it 
to the nullclines of u connected by instanta¬ 
neous horizontal jumps: an approximation that 
becomes exact in the limit of c —> oo. Fur¬ 
ther, to make the model analytically solvable 
we approximate the nullclines of u by linear 
functions through the replacement: /(it, it) —>■ 
fpia{u, v) = v — $(«). Using the instantaneous- 
jump approximation, we obtain v on v = <I>(?x) 
by solving the first order linear equation v = 
1 — v + <b~ 1 (u). 


<I>(it)= 


(Umax — 2oi) 

a—(3 


Um ax 2/3 


p< u< a 
2(X <C U <C Umax 


(18) 


Where we have defined /? = \{u m ax — 

V u max-^/c) and a = We solve 

the equation in the intervals /3 < u < a and 
2 a < u < Umax individually and then switch¬ 
ing between the solutions in the two intervals 


when the solutions reach (it, v) = (a, 2a) and 
(it, u) = (it ma3: , u rnax '). 


As can be seen from Fig. [TJ the resulting 
nullclines of it are linear functions and we are 
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FIG. 8: Left - Numeric solution to the single 
Brusselator at c=50. Right - Solution from the 
PLA for the same parameters. 


now interested in solving: 

= f P ia(u, v) = v - <F(u) (19) 

v = g(u , v) = 1 — v + u 

In order to examine the dynamics close to 
the nullclines, it is more convenient to work in 
terms of u rather than v: ^ ^ = 

_ (ttm a-ff 2a) h! = for 0 < u < a and 

hi = hi for 2a < u < u max- The resulting 
differential equation for u is: 


u = 


— ^[1 + u — $(u)] f3 < u < a 

1 2(X <C U Umax 


( 20 ) 


In the large c regime, 

(2 + S + c) 2 c 

Umax - ^ ~ 4 ’ 

and we know that a > f3 and Umax > 2 a. In 
addition to the time evolution described above, 
the complete dynamics includes instantaneous 
jumps between the two domains: L = (/?, a) 
R = (2a, Umax)- As seen from Fig. |8j the PLA 
captures the functional form of u(t) except for 
the concavity of L. 


agram of coupled oscillators, and thus provide 
additional insight into the mechanisms underly¬ 
ing multistability and a complex attractor land¬ 
scape. In this section, we explicitly study the 
system of two diffusively coupled Brusselators 
using the PLA. The dynamics close to the null¬ 
clines is now generalized to: 


B. PLA analysis of two coupled 
Brusselators 

The real advantage of the PLA is that it al¬ 
lows us to explore the attractors and phase di- 


u'i 


«2 


-i[l - $(ui) + Ml - d(u 1 
1 - d(ui - u 2 ) 

-i[l - $(u 2 ) +u 2 + d(u\ 
1 + d(ui - u 2 ) 


u 2 )} u\ £ L 
ui £ R 

u 2 )] u 2 £ L 
u 2 £ R 


( 21 ) 


( 22 ) 
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The salient feature that emerges from these 
equations is that a Brusselator on the L branch 
experiences an effective coupling that is opposite 
in sign to the coupling it experiences on branch 
R: the coupling is repulsive on L and attractive 
on R. 

In the PLA, the original nonlinearity of the 
equations emerges as jump conditions between 
two distinct branches of linear dynamics. An¬ 
other consequence of this partitioning of the dy¬ 
namics is that the coupling between oscillators 
changes sign between branches. This is the ba¬ 
sic mechanism underlying the complexity in all 
N r rings. 

Since PLA provides a very good approxima¬ 
tion to the dynamics at large c, we use it to 
construct a phase diagram in S — d space by 
defining boundaries between the states listed in 
Table [IJ and comparing it to the one obtained 
from numerical results. The finally comparison 
can be seen in Fig. [lO]We first investigate the 
boundary between amplitude death and out of 
phase oscillations. We note from linear stability 
analysis that for <5 < 0 the fixed point is stable 
and therefore amplitude death can occur. Next 
we examine the bound on when out of phase 
oscillations can occur. To do this we imagine 
both Brusselators to be on the L branch: one 
at (iti = /3) and one just about to jump from 
U 2 = a to «2 = 2 a. 

Ui= -[1 + /3 - Umax + d(a - /3)] (23) 

r? 

ii2 = -[1 — a — d(a — /3)] (24) 

V 

If ii 2 < 0 then the Brusselators will always 
fall to the fixed point. If on the other hand 
ft 2 > 0 then there exist stable out of phase 
oscillations. The phase boundary between be¬ 
tween amplitude death, and multistability of 
amplitude death and out of phase oscillations 
is, therefore, determined by the condition u ,2 = 
0. It is then apparent from Eq. [24] that for 
d > the coupling can cause the pair of 

Brusselators to oscillate. We now examine the 
conditions for the occurrence of in-phase oscil¬ 
lations. We assume the oscillators start out in 


(a) 



FIG. 9: (Color Online) A schematic to 
illustrate the criteria used for constructing the 
phase diagram from PLA: (a) boundary 
between amplitude death and out of phase 
oscillations (b) condition for in-phase 
oscillations, and (c) condition for all 
oscillations being in phase. 


phase and we examine the behavior around the 
u = a to u = 2a jump. We assume the oscil¬ 
lators progress along the nullcline with u\ ~ U 2 
but one oscillator is infinitesimally ahead of the 
other and jumps from a to 2 a first. Without 
loss of generality we say that U 2 jumps before 
«i. Then we have the oscillators at the coordi¬ 
nates (ui, U 2 ) = (a, 2a). At these coordinates 
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FIG. 10: (Color Online) Left - Phase diagram 
constructed from numerical simulations. Right 
- Phase diagram constructed from PLA. 
Significant differences include the sizes of 2 
and 8. Regions 2-7 are d — 6 intervals where 
we expect out of phase oscillations. The left 
bound on this region obtained from PLA is a 
lower bound to the numerical results, since we 
assumed that one of the Brusselators would 
start at a instead of a — d. The PLA, 
therefore, underestimates the distance to be 
travelled along the nullcline, leading to an 
underestimation of the coupling strength d. 
Oscillator death occurs in regions 3, 5, and 7. 
The PLA estimation of these regions should be 
exact. Regions 6-8 are where in-phase 
oscillations can occur. The full nonlinear 
equations increase the stability of in-phase 
oscillations since for a finite c the jump form a 
to 2 a is not instantaneous. 


the derivatives are 

111 = -(1 — a + da) (25) 

V 

11 2 =1 — da (26) 

Now we can see that if d < 1 — — , iii is positive 
implying that the coupling is not strong enough 
to keep ui in L, and therefore, the Brusselators 
can oscillate in phase. 

To find a good bound on where all oscillations 
must be in phase we start with the Brusselators 
out of phase and see when the coupling leads 
to in phase behavior. We first look at the sin¬ 
gle Brusselator dynamics. The assumption here 
is that time scales are not affected stongly by 
the coupling in the regime where the coupling 
is small enough that oscillations are in phase. 
We calculate the time it takes to go along the 
branch L. This is done by solving the first or¬ 
der differential equation with initial conditions 
at u( 0) = j3 and we find that the Brusselator is 
ready to jump at u(t\) = a where to first order 


11 = log W+ij (27) 

The velocity on the right branch is u = 
1 so if the 2 Brusselators start (ui = 
(3,U2 = a) they will progress to ap¬ 

proximately « a, U 2 ~ 2a + log ( 2 ( 1 + 5 )))- 
With derivatives 


iii 


ii2 


=-( 1 — a + d 

V 


=1 — d 


a + log 


a + log 
c 2 > 

2(1 + S)J\ 


2(1 + 6 ) 


(28) 


Now we check when the coupling is strong 
enough to keep Brusselator 1 from jumping 
at these coordinates. We find that for d < 


- a , 1 2 —y the coupling is not strong enough 

«+ lo s(2(fq+J 

to keep both Brusselators off the right branch, 
even if they start out of phase. When they are 
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on the right branch the attractive coupling can 
then pull Brusselators together to form in phase 
oscillations. 

When we combine the above criteria with 
condition for oscillator death discussed in sec¬ 
tion El a phase portrait emerges, which is 
remarkably similar to the one obtained from nu¬ 
merical simulations, as seen from Fig. [TO] 

By greatly simplifying the dynamics we have 
gained some insight into the behavior of the full 
system when the nonlinearity is strong. 

VI. CONCLUSION 

We have presented a detailed analysis of pat¬ 
tern formation in coupled Brusselators in the 
regime of high nonlinearity. This regime is char¬ 
acterized by fast inhibitor dynamics. Coupling 
via this fast species leads to preference for out- 
of-phase oscillations: a feature that is absent 
in the regime of low or moderate nonlinearity. 
Numerical studies of rings of N r coupled oscil¬ 
lators lead to a rich phase diagram with regions 
of oscillator death, in and out of phase oscilla¬ 
tions, and multistability. Analysis of two cou¬ 
pled Brusselators using a piecewise linear ap¬ 
proximation provides a detailed picture of the 
dynamics that leads to out-of-phase oscillations. 
The crucial observation is that in the limit of 
very strong nonlinearity the limit cycle can be 
broken down into two branches in one of which 


the coupling is repulsive and in the other the 
coupling is attractive. A phase diagram based 
on the PLA is found to be in good agreement 
with numerical results. 

Numerical analysis for a ring of 5 oscillators 
yields patterns in which the frequencies of the 
oscillators synchronize with integer frequency 
ratios. For rings with fewer than 5 oscilla¬ 
tors, the diversity of patterns is characterized by 
their phase relationships since the frequencies of 
all oscillators are identical. The PLA analysis 
suggests that one can reduce the coupled dy¬ 
namics of strongly nonlinear Brusselators to a 
discrete-time map, which are know to exhibit 
frequency lockmg[21]. This analysis and a de¬ 
tailed numerical study of the complete phase 
diagram for the N r = 3, and 4 rings will be the 
subject of future research. 
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